﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.Arm;
using System.Runtime.Intrinsics.X86;
using System.Text;
using System.Threading.Tasks;

namespace VectorCalculator
{
    internal class NormalizationWithSIMD
    {

        private const int AVX_VECTOR_SIZE = 8;   // AVX2处理8个float
        private const int NEON_VECTOR_SIZE = 4;  // NEON处理4个float
        private const int SSE_VECTOR_SIZE = 4;   // SSE处理4个float


        #region 向量归一化（支持多算法 + SIMD）

        public static float[] Normalize(float[] vector, NormalizationType type = NormalizationType.L2)
        {
            if (vector == null || vector.Length == 0)
                throw new ArgumentException("向量不能为空");

            var result = new float[vector.Length];

            switch (type)
            {
                case NormalizationType.L2:
                    Normalize_L2(vector, result);
                    break;
                case NormalizationType.L1:
                    Normalize_L1(vector, result);
                    break;
                case NormalizationType.Max:
                    Normalize_Max(vector, result);
                    break;
                case NormalizationType.MinMax:
                    Normalize_MinMax(vector, result);
                    break;
                case NormalizationType.ZScore:
                    Normalize_ZScore(vector, result);
                    break;
                default:
                    throw new NotSupportedException($"不支持的归一化类型：{type}");
            }

            return result;
        }

        #region L2

        private static void Normalize_L2(float[] vector, float[] result)
        {
            float norm = ComputeL2Norm(vector);
            if (norm < float.Epsilon)
            {
                // 零向量：按常见做法直接返回原向量（也可返回0向量，按需改）
                Array.Copy(vector, result, vector.Length);
                return;
            }
            ScaleWrite(vector, result, 1.0f / norm);
        }

        #endregion

        #region L1（绝对值和为1）

        private static void Normalize_L1(float[] vector, float[] result)
        {
            float sumAbs = ComputeSumAbs(vector);
            if (sumAbs < float.Epsilon)
            {
                Array.Copy(vector, result, vector.Length);
                return;
            }
            ScaleWrite(vector, result, 1.0f / sumAbs);
        }

        #endregion

        #region Max（最大绝对值为1）

        private static void Normalize_Max(float[] vector, float[] result)
        {
            float maxAbs = ComputeMaxAbs(vector);
            if (maxAbs < float.Epsilon)
            {
                Array.Copy(vector, result, vector.Length);
                return;
            }
            ScaleWrite(vector, result, 1.0f / maxAbs);
        }

        #endregion

        #region Min-Max（映射到[0,1]）

        private static void Normalize_MinMax(float[] vector, float[] result)
        {
            (float min, float max) = ComputeMinMax(vector);
            float range = max - min;
            if (range < float.Epsilon)
            {
                // 常数向量，全部置0（或置1/0.5，按业务需要改）
                Array.Clear(result, 0, result.Length);
                return;
            }

            if (Avx.IsSupported && vector.Length >= AVX_VECTOR_SIZE)
                AffineWriteAvx(vector, result, -min, 1.0f / range); // (x - min) / range == x*(1/range) + (-min)*(1/range)
            else if (AdvSimd.IsSupported && vector.Length >= NEON_VECTOR_SIZE)
                AffineWriteNeon(vector, result, -min, 1.0f / range);
            else if (Sse.IsSupported && vector.Length >= SSE_VECTOR_SIZE)
                AffineWriteSse(vector, result, -min, 1.0f / range);
            else
                AffineWriteScalar(vector, result, -min, 1.0f / range);
        }

        #endregion

        #region Z-Score（均值0，标准差1）

        private static void Normalize_ZScore(float[] vector, float[] result)
        {
            (double mean, double std) = ComputeMeanStd(vector);
            if (std < 1e-12)
            {
                // 标准差接近0：返回全0（或全1），这里返回0
                Array.Clear(result, 0, result.Length);
                return;
            }

            float fMean = (float)mean;
            float invStd = (float)(1.0 / std);

            if (Avx.IsSupported && vector.Length >= AVX_VECTOR_SIZE)
                ZScoreWriteAvx(vector, result, fMean, invStd);
            else if (AdvSimd.IsSupported && vector.Length >= NEON_VECTOR_SIZE)
                ZScoreWriteNeon(vector, result, fMean, invStd);
            else if (Sse.IsSupported && vector.Length >= SSE_VECTOR_SIZE)
                ZScoreWriteSse(vector, result, fMean, invStd);
            else
                ZScoreWriteScalar(vector, result, fMean, invStd);
        }

        #endregion

        /*==================== 统计阶段（SIMD） ====================*/

        private static unsafe float ComputeL2Norm(float[] vector)
        {
            float sumSq = 0f;

            if (Avx.IsSupported && vector.Length >= AVX_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % AVX_VECTOR_SIZE);
                Vector256<float> acc = Vector256<float>.Zero;
                fixed (float* p = vector)
                {
                    for (int i = 0; i < m; i += AVX_VECTOR_SIZE)
                    {
                        var v = Avx.LoadVector256(p + i);
                        acc = Avx.Add(acc, Avx.Multiply(v, v));
                    }
                }
                sumSq = HorizontalSum(acc);
                for (int i = m; i < vector.Length; i++) sumSq += vector[i] * vector[i];
                return MathF.Sqrt(sumSq);
            }
            else if (AdvSimd.IsSupported && vector.Length >= NEON_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % NEON_VECTOR_SIZE);
                Vector128<float> acc = Vector128<float>.Zero;
                fixed (float* p = vector)
                {
                    for (int i = 0; i < m; i += NEON_VECTOR_SIZE)
                    {
                        var v = AdvSimd.LoadVector128(p + i);
                        acc = AdvSimd.Add(acc, AdvSimd.Multiply(v, v));
                    }
                }
                sumSq = HorizontalSum(acc);
                for (int i = m; i < vector.Length; i++) sumSq += vector[i] * vector[i];
                return MathF.Sqrt(sumSq);
            }
            else if (Sse.IsSupported && vector.Length >= SSE_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % SSE_VECTOR_SIZE);
                Vector128<float> acc = Vector128<float>.Zero;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        for (int i = 0; i < m; i += SSE_VECTOR_SIZE)
                        {
                            var v = Sse.LoadVector128(p + i);
                            acc = Sse.Add(acc, Sse.Multiply(v, v));
                        }
                    }
                }
                sumSq = HorizontalSum(acc);
                for (int i = m; i < vector.Length; i++) sumSq += vector[i] * vector[i];
                return MathF.Sqrt(sumSq);
            }
            else
            {
                for (int i = 0; i < vector.Length; i++) sumSq += vector[i] * vector[i];
                return MathF.Sqrt(sumSq);
            }
        }

        private static unsafe float ComputeSumAbs(float[] vector)
        {
            float sum = 0f;
            if (Avx.IsSupported && vector.Length >= AVX_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % AVX_VECTOR_SIZE);
                var signMask = Vector256.Create(-0.0f);
                Vector256<float> acc = Vector256<float>.Zero;
                fixed (float* p = vector)
                {
                    for (int i = 0; i < m; i += AVX_VECTOR_SIZE)
                    {
                        var v = Avx.LoadVector256(p + i);
                        var abs = Avx.AndNot(signMask, v); // 清符号位
                        acc = Avx.Add(acc, abs);
                    }
                }
                sum = HorizontalSum(acc);
                for (int i = m; i < vector.Length; i++) sum += MathF.Abs(vector[i]);
                return sum;
            }
            else if (AdvSimd.IsSupported && vector.Length >= NEON_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % NEON_VECTOR_SIZE);
                Vector128<float> acc = Vector128<float>.Zero;
                fixed (float* p = vector)
                {
                    for (int i = 0; i < m; i += NEON_VECTOR_SIZE)
                    {
                        var v = AdvSimd.LoadVector128(p + i);
                        var abs = AdvSimd.Abs(v);
                        acc = AdvSimd.Add(acc, abs);
                    }
                }
                sum = HorizontalSum(acc);
                for (int i = m; i < vector.Length; i++) sum += MathF.Abs(vector[i]);
                return sum;
            }
            else if (Sse.IsSupported && vector.Length >= SSE_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % SSE_VECTOR_SIZE);
                var signMask = Vector128.Create(-0.0f);
                Vector128<float> acc = Vector128<float>.Zero;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        for (int i = 0; i < m; i += SSE_VECTOR_SIZE)
                        {
                            var v = Sse.LoadVector128(p + i);
                            var abs = Sse.AndNot(signMask, v);
                            acc = Sse.Add(acc, abs);
                        }
                    }
                }
                sum = HorizontalSum(acc);
                for (int i = m; i < vector.Length; i++) sum += MathF.Abs(vector[i]);
                return sum;
            }
            else
            {
                for (int i = 0; i < vector.Length; i++) sum += MathF.Abs(vector[i]);
                return sum;
            }
        }

        private static unsafe float ComputeMaxAbs(float[] vector)
        {
            float maxAbs = 0f;
            if (Avx.IsSupported && vector.Length >= AVX_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % AVX_VECTOR_SIZE);
                var signMask = Vector256.Create(-0.0f);
                Vector256<float> acc = Vector256<float>.Zero; // 先用0
                bool accInit = false;

                fixed (float* p = vector)
                {
                    for (int i = 0; i < m; i += AVX_VECTOR_SIZE)
                    {
                        var v = Avx.LoadVector256(p + i);
                        var abs = Avx.AndNot(signMask, v);
                        if (!accInit) { acc = abs; accInit = true; }
                        else { acc = Avx.Max(acc, abs); }
                    }
                }
                maxAbs = HorizontalMax(acc);
                for (int i = m; i < vector.Length; i++) maxAbs = MathF.Max(maxAbs, MathF.Abs(vector[i]));
                return maxAbs;
            }
            else if (AdvSimd.IsSupported && vector.Length >= NEON_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % NEON_VECTOR_SIZE);
                Vector128<float> acc = Vector128<float>.Zero;
                bool accInit = false;
                fixed (float* p = vector)
                {
                    for (int i = 0; i < m; i += NEON_VECTOR_SIZE)
                    {
                        var v = AdvSimd.LoadVector128(p + i);
                        var abs = AdvSimd.Abs(v);
                        acc = accInit ? AdvSimd.Max(acc, abs) : abs;
                        accInit = true;
                    }
                }
                maxAbs = HorizontalMax(acc);
                for (int i = m; i < vector.Length; i++) maxAbs = MathF.Max(maxAbs, MathF.Abs(vector[i]));
                return maxAbs;
            }
            else if (Sse.IsSupported && vector.Length >= SSE_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % SSE_VECTOR_SIZE);
                var signMask = Vector128.Create(-0.0f);
                Vector128<float> acc = Vector128<float>.Zero;
                bool accInit = false;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        for (int i = 0; i < m; i += SSE_VECTOR_SIZE)
                        {
                            var v = Sse.LoadVector128(p + i);
                            var abs = Sse.AndNot(signMask, v);
                            acc = accInit ? Sse.Max(acc, abs) : abs;
                            accInit = true;
                        }
                    }
                }
                maxAbs = HorizontalMax(acc);
                for (int i = m; i < vector.Length; i++) maxAbs = MathF.Max(maxAbs, MathF.Abs(vector[i]));
                return maxAbs;
            }
            else
            {
                for (int i = 0; i < vector.Length; i++) maxAbs = MathF.Max(maxAbs, MathF.Abs(vector[i]));
                return maxAbs;
            }
        }

        private static unsafe (float min, float max) ComputeMinMax(float[] vector)
        {
            if (vector.Length == 0) return (0, 0);

            float min = vector[0], max = vector[0];

            if (Avx.IsSupported && vector.Length >= AVX_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % AVX_VECTOR_SIZE);
                Vector256<float> vMin, vMax;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        vMin = Avx.LoadVector256(p);
                        vMax = vMin;
                        for (int i = AVX_VECTOR_SIZE; i < m; i += AVX_VECTOR_SIZE)
                        {
                            var v = Avx.LoadVector256(p + i);
                            vMin = Avx.Min(vMin, v);
                            vMax = Avx.Max(vMax, v);
                        }
                    }
                }
                min = HorizontalMin(vMin);
                max = HorizontalMax(vMax);
                for (int i = m; i < vector.Length; i++)
                {
                    min = MathF.Min(min, vector[i]);
                    max = MathF.Max(max, vector[i]);
                }
                return (min, max);
            }
            else if (AdvSimd.IsSupported && vector.Length >= NEON_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % NEON_VECTOR_SIZE);
                Vector128<float> vMin, vMax;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        vMin = AdvSimd.LoadVector128(p);
                        vMax = vMin;
                        for (int i = NEON_VECTOR_SIZE; i < m; i += NEON_VECTOR_SIZE)
                        {
                            var v = AdvSimd.LoadVector128(p + i);
                            vMin = AdvSimd.Min(vMin, v);
                            vMax = AdvSimd.Max(vMax, v);
                        }
                    }
                }
                min = HorizontalMin(vMin);
                max = HorizontalMax(vMax);
                for (int i = m; i < vector.Length; i++)
                {
                    min = MathF.Min(min, vector[i]);
                    max = MathF.Max(max, vector[i]);
                }
                return (min, max);
            }
            else if (Sse.IsSupported && vector.Length >= SSE_VECTOR_SIZE)
            {
                int n = vector.Length, m = n - (n % SSE_VECTOR_SIZE);
                Vector128<float> vMin, vMax;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        vMin = Sse.LoadVector128(p);
                        vMax = vMin;
                        for (int i = SSE_VECTOR_SIZE; i < m; i += SSE_VECTOR_SIZE)
                        {
                            var v = Sse.LoadVector128(p + i);
                            vMin = Sse.Min(vMin, v);
                            vMax = Sse.Max(vMax, v);
                        }
                    }
                }
                min = HorizontalMin(vMin);
                max = HorizontalMax(vMax);
                for (int i = m; i < vector.Length; i++)
                {
                    min = MathF.Min(min, vector[i]);
                    max = MathF.Max(max, vector[i]);
                }
                return (min, max);
            }
            else
            {
                for (int i = 1; i < vector.Length; i++)
                {
                    float v = vector[i];
                    if (v < min) min = v;
                    if (v > max) max = v;
                }
                return (min, max);
            }
        }

        private static unsafe (double mean, double std) ComputeMeanStd(float[] vector)
        {
            // 同时累加 sum 和 sumsq，减少一次内存遍历
            double sum = 0.0, sumSq = 0.0;
            int n = vector.Length;

            if (Avx.IsSupported && n >= AVX_VECTOR_SIZE)
            {
                int m = n - (n % AVX_VECTOR_SIZE);
                Vector256<float> acc1 = Vector256<float>.Zero;
                Vector256<float> acc2 = Vector256<float>.Zero;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        for (int i = 0; i < m; i += AVX_VECTOR_SIZE)
                        {
                            var v = Avx.LoadVector256(p + i);
                            acc1 = Avx.Add(acc1, v);
                            acc2 = Avx.Add(acc2, Avx.Multiply(v, v));
                        }
                    }
                }
                sum += HorizontalSum(acc1);
                sumSq += HorizontalSum(acc2);
                for (int i = m; i < n; i++)
                {
                    float v = vector[i];
                    sum += v;
                    sumSq += v * v;
                }
            }
            else if (AdvSimd.IsSupported && n >= NEON_VECTOR_SIZE)
            {
                int m = n - (n % NEON_VECTOR_SIZE);
                Vector128<float> acc1 = Vector128<float>.Zero;
                Vector128<float> acc2 = Vector128<float>.Zero;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        for (int i = 0; i < m; i += NEON_VECTOR_SIZE)
                        {
                            var v = AdvSimd.LoadVector128(p + i);
                            acc1 = AdvSimd.Add(acc1, v);
                            acc2 = AdvSimd.Add(acc2, AdvSimd.Multiply(v, v));
                        }
                    }
                }
                sum += HorizontalSum(acc1);
                sumSq += HorizontalSum(acc2);
                for (int i = m; i < n; i++)
                {
                    float v = vector[i];
                    sum += v;
                    sumSq += v * v;
                }
            }
            else if (Sse.IsSupported && n >= SSE_VECTOR_SIZE)
            {
                int m = n - (n % SSE_VECTOR_SIZE);
                Vector128<float> acc1 = Vector128<float>.Zero;
                Vector128<float> acc2 = Vector128<float>.Zero;
                unsafe
                {
                    fixed (float* p = vector)
                    {
                        for (int i = 0; i < m; i += SSE_VECTOR_SIZE)
                        {
                            var v = Sse.LoadVector128(p + i);
                            acc1 = Sse.Add(acc1, v);
                            acc2 = Sse.Add(acc2, Sse.Multiply(v, v));
                        }
                    }
                }
                sum += HorizontalSum(acc1);
                sumSq += HorizontalSum(acc2);
                for (int i = m; i < n; i++)
                {
                    float v = vector[i];
                    sum += v;
                    sumSq += v * v;
                }
            }
            else
            {
                for (int i = 0; i < n; i++)
                {
                    float v = vector[i];
                    sum += v;
                    sumSq += v * v;
                }
            }

            double mean = sum / Math.Max(1, n);
            // 方差的数值稳定性：使用 E[x^2] - (E[x])^2
            double var = Math.Max(0.0, (sumSq / Math.Max(1, n)) - mean * mean);
            double std = Math.Sqrt(var);
            return (mean, std);
        }

        /*==================== 写回阶段（SIMD） ====================*/

        private static void ScaleWrite(float[] vector, float[] result, float scale)
        {
            if (Avx.IsSupported && vector.Length >= AVX_VECTOR_SIZE)
                ScaleWriteAvx(vector, result, scale);
            else if (AdvSimd.IsSupported && vector.Length >= NEON_VECTOR_SIZE)
                ScaleWriteNeon(vector, result, scale);
            else if (Sse.IsSupported && vector.Length >= SSE_VECTOR_SIZE)
                ScaleWriteSse(vector, result, scale);
            else
                ScaleWriteScalar(vector, result, scale);
        }

        private static unsafe void ScaleWriteAvx(float[] vector, float[] result, float scale)
        {
            int n = vector.Length, m = n - (n % AVX_VECTOR_SIZE);
            var s = Vector256.Create(scale);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += AVX_VECTOR_SIZE)
                {
                    var v = Avx.LoadVector256(pv + i);
                    var w = Avx.Multiply(v, s);
                    Avx.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = vector[i] * scale;
        }

        private static unsafe void ScaleWriteSse(float[] vector, float[] result, float scale)
        {
            int n = vector.Length, m = n - (n % SSE_VECTOR_SIZE);
            var s = Vector128.Create(scale);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += SSE_VECTOR_SIZE)
                {
                    var v = Sse.LoadVector128(pv + i);
                    var w = Sse.Multiply(v, s);
                    Sse.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = vector[i] * scale;
        }

        private static unsafe void ScaleWriteNeon(float[] vector, float[] result, float scale)
        {
            int n = vector.Length, m = n - (n % NEON_VECTOR_SIZE);
            var s = Vector128.Create(scale);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += NEON_VECTOR_SIZE)
                {
                    var v = AdvSimd.LoadVector128(pv + i);
                    var w = AdvSimd.Multiply(v, s);
                    AdvSimd.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = vector[i] * scale;
        }

        private static void ScaleWriteScalar(float[] vector, float[] result, float scale)
        {
            for (int i = 0; i < vector.Length; i++) result[i] = vector[i] * scale;
        }

        // y = (x + add) * mul
        private static void AffineWriteScalar(float[] vector, float[] result, float add, float mul)
        {
            for (int i = 0; i < vector.Length; i++)
                result[i] = (vector[i] + add) * mul;
        }

        private static unsafe void AffineWriteAvx(float[] vector, float[] result, float add, float mul)
        {
            int n = vector.Length, m = n - (n % AVX_VECTOR_SIZE);
            var vAdd = Vector256.Create(add);
            var vMul = Vector256.Create(mul);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += AVX_VECTOR_SIZE)
                {
                    var v = Avx.LoadVector256(pv + i);
                    var t = Avx.Add(v, vAdd);
                    var w = Avx.Multiply(t, vMul);
                    Avx.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = (vector[i] + add) * mul;
        }

        private static unsafe void AffineWriteSse(float[] vector, float[] result, float add, float mul)
        {
            int n = vector.Length, m = n - (n % SSE_VECTOR_SIZE);
            var vAdd = Vector128.Create(add);
            var vMul = Vector128.Create(mul);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += SSE_VECTOR_SIZE)
                {
                    var v = Sse.LoadVector128(pv + i);
                    var t = Sse.Add(v, vAdd);
                    var w = Sse.Multiply(t, vMul);
                    Sse.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = (vector[i] + add) * mul;
        }

        private static unsafe void AffineWriteNeon(float[] vector, float[] result, float add, float mul)
        {
            int n = vector.Length, m = n - (n % NEON_VECTOR_SIZE);
            var vAdd = Vector128.Create(add);
            var vMul = Vector128.Create(mul);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += NEON_VECTOR_SIZE)
                {
                    var v = AdvSimd.LoadVector128(pv + i);
                    var t = AdvSimd.Add(v, vAdd);
                    var w = AdvSimd.Multiply(t, vMul);
                    AdvSimd.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = (vector[i] + add) * mul;
        }

        // y = (x - mean) * invStd
        private static void ZScoreWriteScalar(float[] vector, float[] result, float mean, float invStd)
        {
            for (int i = 0; i < vector.Length; i++)
                result[i] = (vector[i] - mean) * invStd;
        }

        private static unsafe void ZScoreWriteAvx(float[] vector, float[] result, float mean, float invStd)
        {
            int n = vector.Length, m = n - (n % AVX_VECTOR_SIZE);
            var vMean = Vector256.Create(mean);
            var vInv = Vector256.Create(invStd);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += AVX_VECTOR_SIZE)
                {
                    var v = Avx.LoadVector256(pv + i);
                    var t = Avx.Subtract(v, vMean);
                    var w = Avx.Multiply(t, vInv);
                    Avx.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = (vector[i] - mean) * invStd;
        }

        private static unsafe void ZScoreWriteSse(float[] vector, float[] result, float mean, float invStd)
        {
            int n = vector.Length, m = n - (n % SSE_VECTOR_SIZE);
            var vMean = Vector128.Create(mean);
            var vInv = Vector128.Create(invStd);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += SSE_VECTOR_SIZE)
                {
                    var v = Sse.LoadVector128(pv + i);
                    var t = Sse.Subtract(v, vMean);
                    var w = Sse.Multiply(t, vInv);
                    Sse.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = (vector[i] - mean) * invStd;
        }

        private static unsafe void ZScoreWriteNeon(float[] vector, float[] result, float mean, float invStd)
        {
            int n = vector.Length, m = n - (n % NEON_VECTOR_SIZE);
            var vMean = Vector128.Create(mean);
            var vInv = Vector128.Create(invStd);
            fixed (float* pv = vector, pr = result)
            {
                for (int i = 0; i < m; i += NEON_VECTOR_SIZE)
                {
                    var v = AdvSimd.LoadVector128(pv + i);
                    var t = AdvSimd.Subtract(v, vMean);
                    var w = AdvSimd.Multiply(t, vInv);
                    AdvSimd.Store(pr + i, w);
                }
            }
            for (int i = m; i < n; i++) result[i] = (vector[i] - mean) * invStd;
        }

        /*==================== 水平归约辅助 ====================*/
        // 方法内部使用 unsafe 块

        private static float HorizontalSum(Vector256<float> v)
        {
            unsafe
            {
                // 存到栈上然后逐元素相加
                float* tmp = stackalloc float[8];
                Avx.Store(tmp, v);
                return tmp[0] + tmp[1] + tmp[2] + tmp[3] + tmp[4] + tmp[5] + tmp[6] + tmp[7];
            }
        }

        private static float HorizontalSum(Vector128<float> v)
        {
            unsafe
            {
                float* tmp = stackalloc float[4];
                Sse.Store(tmp, v);
                return tmp[0] + tmp[1] + tmp[2] + tmp[3];
            }
        }

        private static float HorizontalMax(Vector256<float> v)
        {
            unsafe
            {
                float* tmp = stackalloc float[8];
                Avx.Store(tmp, v);
                float m = tmp[0];
                // 展开比较，比用循环更快且更明确
                if (tmp[1] > m) m = tmp[1];
                if (tmp[2] > m) m = tmp[2];
                if (tmp[3] > m) m = tmp[3];
                if (tmp[4] > m) m = tmp[4];
                if (tmp[5] > m) m = tmp[5];
                if (tmp[6] > m) m = tmp[6];
                if (tmp[7] > m) m = tmp[7];
                return m;
            }
        }

        private static float HorizontalMax(Vector128<float> v)
        {
            unsafe
            {
                float* tmp = stackalloc float[4];
                Sse.Store(tmp, v);
                float m = tmp[0];
                if (tmp[1] > m) m = tmp[1];
                if (tmp[2] > m) m = tmp[2];
                if (tmp[3] > m) m = tmp[3];
                return m;
            }
        }

        private static float HorizontalMin(Vector256<float> v)
        {
            unsafe
            {
                float* tmp = stackalloc float[8];
                Avx.Store(tmp, v);
                float m = tmp[0];
                if (tmp[1] < m) m = tmp[1];
                if (tmp[2] < m) m = tmp[2];
                if (tmp[3] < m) m = tmp[3];
                if (tmp[4] < m) m = tmp[4];
                if (tmp[5] < m) m = tmp[5];
                if (tmp[6] < m) m = tmp[6];
                if (tmp[7] < m) m = tmp[7];
                return m;
            }
        }

        private static float HorizontalMin(Vector128<float> v)
        {
            unsafe
            {
                float* tmp = stackalloc float[4];
                Sse.Store(tmp, v);
                float m = tmp[0];
                if (tmp[1] < m) m = tmp[1];
                if (tmp[2] < m) m = tmp[2];
                if (tmp[3] < m) m = tmp[3];
                return m;
            }
        }


        #endregion

    }
}
